Usenet 1993 July

home *** CD-ROM | disk | FTP | other *** search

/ Usenet 1993 July / InfoMagic USENET CD-ROM July 1993.ISO / sources / misc / volume24 / gnucalc / part05 < prev next >

Wrap

Text File | 1991-10-28 | 55.5 KB | 1,782 lines

Newsgroups: comp.sources.misc From: daveg@synaptics.com (David Gillespie) Subject: v24i053: gnucalc - GNU Emacs Calculator, v2.00, Part05/56 Message-ID: <1991Oct29.042451.7188@sparky.imd.sterling.com> X-Md4-Signature: c9449dd13a7867353f1d8e87767c2a9c Date: Tue, 29 Oct 1991 04:24:51 GMT Approved: kent@sparky.imd.sterling.com Submitted-by: daveg@synaptics.com (David Gillespie) Posting-number: Volume 24, Issue 53 Archive-name: gnucalc/part05 Environment: Emacs Supersedes: gmcalc: Volume 13, Issue 27-45 ---- Cut Here and unpack ---- #!/bin/sh # this is Part.05 (part 5 of a multipart archive) # do not concatenate these parts, unpack them in order with /bin/sh # file calc-alg-2.el continued # if test ! -r _shar_seq_.tmp; then echo 'Please unpack part 1 first!' exit 1 fi (read Scheck if test "$Scheck" != 5; then echo Please unpack part "$Scheck" next! exit 1 else exit 0 fi ) < _shar_seq_.tmp || exit 1 if test ! -f _shar_wnt_.tmp; then echo 'x - still skipping calc-alg-2.el' else echo 'x - continuing file calc-alg-2.el' sed 's/^X//' << 'SHAR_EOF' >> 'calc-alg-2.el' && X (not (math-expr-contains (nth 2 expr) var))) X (math-mul (calcFunc-integ (nth 1 expr) var low high) (nth 2 expr))) X ((and (eq (car-safe expr) '/) X (not (math-expr-contains (nth 1 expr) var)) X (not (math-equal-int (nth 1 expr) 1))) X (math-mul (nth 1 expr) X (calcFunc-integ (math-div 1 (nth 2 expr)) var low high))) X ((and (eq (car-safe expr) '/) X (not (math-expr-contains (nth 2 expr) var))) X (math-div (calcFunc-integ (nth 1 expr) var low high) (nth 2 expr))) X ((eq (car-safe expr) 'vec) X (cons 'vec (mapcar (function (lambda (x) (calcFunc-integ x var low high))) X (cdr expr)))) X (t X (let ((state (list calc-angle-mode X ;;calc-symbolic-mode X ;;calc-prefer-frac X calc-internal-prec X (calc-var-value 'var-IntegRules) X (calc-var-value 'var-IntegSimpRules)))) X (or (equal state math-integral-cache-state) X (setq math-integral-cache-state state X math-integral-cache nil))) X (let* ((math-max-integral-limit (or (and (boundp 'var-IntegLimit) X (natnump var-IntegLimit) X var-IntegLimit) X 3)) X (math-integral-limit 1) X (sexpr (math-expr-subst expr var math-integ-var)) X (trace-buffer (get-buffer "*Trace*")) X (calc-language (if (eq calc-language 'big) nil calc-language)) X (math-any-substs t) X (math-enable-subst nil) X (math-prev-parts-v nil) X (math-doing-parts nil) X (math-good-parts nil) X (res X (if trace-buffer X (let ((calcbuf (current-buffer)) X (calcwin (selected-window))) X (unwind-protect X (progn X (if (get-buffer-window trace-buffer) X (select-window (get-buffer-window trace-buffer))) X (set-buffer trace-buffer) X (goto-char (point-max)) X (or (assq 'scroll-stop (buffer-local-variables)) X (progn X (make-local-variable 'scroll-step) X (setq scroll-step 3))) X (insert "\n\n\n") X (set-buffer calcbuf) X (math-try-integral sexpr)) X (select-window calcwin) X (set-buffer calcbuf))) X (math-try-integral sexpr)))) X (if res X (progn X (if (calc-has-rules 'var-IntegAfterRules) X (setq res (math-rewrite res '(var IntegAfterRules X var-IntegAfterRules)))) X (math-simplify X (if (and low high) X (math-sub (math-expr-subst res math-integ-var high) X (math-expr-subst res math-integ-var low)) X (setq res (math-fix-const-terms res math-integ-vars)) X (if low X (math-expr-subst res math-integ-var low) X (math-expr-subst res math-integ-var var))))) X (append (list 'calcFunc-integ expr var) X (and low (list low)) X (and high (list high))))))) ) X X (math-defintegral calcFunc-inv X (math-integral (math-div 1 u))) X (math-defintegral calcFunc-conj X (let ((int (math-integral u))) X (and int X (list 'calcFunc-conj int)))) X (math-defintegral calcFunc-deg X (let ((int (math-integral u))) X (and int X (list 'calcFunc-deg int)))) X (math-defintegral calcFunc-rad X (let ((int (math-integral u))) X (and int X (list 'calcFunc-rad int)))) X (math-defintegral calcFunc-re X (let ((int (math-integral u))) X (and int X (list 'calcFunc-re int)))) X (math-defintegral calcFunc-im X (let ((int (math-integral u))) X (and int X (list 'calcFunc-im int)))) X (math-defintegral calcFunc-sqrt X (and (equal u math-integ-var) X (math-mul '(frac 2 3) X (list 'calcFunc-sqrt (math-pow u 3))))) X (math-defintegral calcFunc-exp X (and (equal u math-integ-var) X (list 'calcFunc-exp u))) X (math-defintegral calcFunc-ln X (or (and (equal u math-integ-var) X (math-sub (math-mul u (list 'calcFunc-ln u)) u)) X (and (eq (car u) '*) X (math-integral (math-add (list 'calcFunc-ln (nth 1 u)) X (list 'calcFunc-ln (nth 2 u))))) X (and (eq (car u) '/) X (math-integral (math-sub (list 'calcFunc-ln (nth 1 u)) X (list 'calcFunc-ln (nth 2 u))))) X (and (eq (car u) '^) X (math-integral (math-mul (nth 2 u) X (list 'calcFunc-ln (nth 1 u))))))) X (math-defintegral calcFunc-log10 X (and (equal u math-integ-var) X (math-sub (math-mul u (list 'calcFunc-ln u)) X (math-div u (list 'calcFunc-ln 10))))) X (math-defintegral-2 calcFunc-log X (math-integral (math-div (list 'calcFunc-ln u) X (list 'calcFunc-ln v)))) X (math-defintegral calcFunc-sin X (and (equal u math-integ-var) X (math-neg (math-from-radians-2 (list 'calcFunc-cos u))))) X (math-defintegral calcFunc-cos X (and (equal u math-integ-var) X (math-from-radians-2 (list 'calcFunc-sin u)))) X (math-defintegral calcFunc-tan X (and (equal u math-integ-var) X (math-neg (math-from-radians-2 X (list 'calcFunc-ln (list 'calcFunc-cos u)))))) X (math-defintegral calcFunc-arcsin X (and (equal u math-integ-var) X (math-add (math-mul u (list 'calcFunc-arcsin u)) X (math-from-radians-2 X (list 'calcFunc-sqrt (math-sub 1 (math-sqr u))))))) X (math-defintegral calcFunc-arccos X (and (equal u math-integ-var) X (math-sub (math-mul u (list 'calcFunc-arccos u)) X (math-from-radians-2 X (list 'calcFunc-sqrt (math-sub 1 (math-sqr u))))))) X (math-defintegral calcFunc-arctan X (and (equal u math-integ-var) X (math-sub (math-mul u (list 'calcFunc-arctan u)) X (math-from-radians-2 X (math-div (list 'calcFunc-ln (math-add 1 (math-sqr u))) X 2))))) X (math-defintegral calcFunc-sinh X (and (equal u math-integ-var) X (list 'calcFunc-cosh u))) X (math-defintegral calcFunc-cosh X (and (equal u math-integ-var) X (list 'calcFunc-sinh u))) X (math-defintegral calcFunc-tanh X (and (equal u math-integ-var) X (list 'calcFunc-ln (list 'calcFunc-cosh u)))) X (math-defintegral calcFunc-arcsinh X (and (equal u math-integ-var) X (math-sub (math-mul u (list 'calcFunc-arcsinh u)) X (list 'calcFunc-sqrt (math-add (math-sqr u) 1))))) X (math-defintegral calcFunc-arccosh X (and (equal u math-integ-var) X (math-sub (math-mul u (list 'calcFunc-arccosh u)) X (list 'calcFunc-sqrt (math-sub 1 (math-sqr u)))))) X (math-defintegral calcFunc-arctanh X (and (equal u math-integ-var) X (math-sub (math-mul u (list 'calcFunc-arctan u)) X (math-div (list 'calcFunc-ln X (math-add 1 (math-sqr u))) X 2)))) X ;;; (Ax + B) / (ax^2 + bx + c)^n forms. (math-defintegral-2 / X (math-integral-rational-funcs u v)) X (defun math-integral-rational-funcs (u v) X (let ((pu (math-is-polynomial u math-integ-var 1)) X (vpow 1) pv) X (and pu X (catch 'int-rat X (if (and (eq (car-safe v) '^) (natnump (nth 2 v))) X (setq vpow (nth 2 v) X v (nth 1 v))) X (and (setq pv (math-is-polynomial v math-integ-var 2)) X (let ((int (math-mul-thru X (car pu) X (math-integral-q02 (car pv) (nth 1 pv) X (nth 2 pv) v vpow)))) X (if (cdr pu) X (setq int (math-add int X (math-mul-thru X (nth 1 pu) X (math-integral-q12 X (car pv) (nth 1 pv) X (nth 2 pv) v vpow))))) X int)))))) X (defun math-integral-q12 (a b c v vpow) X (let (q) X (cond ((not c) X (cond ((= vpow 1) X (math-sub (math-div math-integ-var b) X (math-mul (math-div a (math-sqr b)) X (list 'calcFunc-ln v)))) X ((= vpow 2) X (math-div (math-add (list 'calcFunc-ln v) X (math-div a v)) X (math-sqr b))) X (t X (let ((nm1 (math-sub vpow 1)) X (nm2 (math-sub vpow 2))) X (math-div (math-sub X (math-div a (math-mul nm1 (math-pow v nm1))) X (math-div 1 (math-mul nm2 (math-pow v nm2)))) X (math-sqr b)))))) X ((math-zerop X (setq q (math-sub (math-mul 4 (math-mul a c)) (math-sqr b)))) X (let ((part (math-div b (math-mul 2 c)))) X (math-mul-thru (math-pow c vpow) X (math-integral-q12 part 1 nil X (math-add math-integ-var part) X (* vpow 2))))) X ((= vpow 1) X (and (math-ratp q) (math-negp q) X (let ((calc-symbolic-mode t)) X (math-ratp (math-sqrt (math-neg q)))) X (throw 'int-rat nil)) ; should have used calcFunc-apart first X (math-sub (math-div (list 'calcFunc-ln v) (math-mul 2 c)) X (math-mul-thru (math-div b (math-mul 2 c)) X (math-integral-q02 a b c v 1)))) X (t X (let ((n (1- vpow))) X (math-sub (math-neg (math-div X (math-add (math-mul b math-integ-var) X (math-mul 2 a)) X (math-mul n (math-mul q (math-pow v n))))) X (math-mul-thru (math-div (math-mul b (1- (* 2 n))) X (math-mul n q)) X (math-integral-q02 a b c v n))))))) ) X (defun math-integral-q02 (a b c v vpow) X (let (q rq part) X (cond ((not c) X (cond ((= vpow 1) X (math-div (list 'calcFunc-ln v) b)) X (t X (math-div (math-pow v (- 1 vpow)) X (math-mul (- 1 vpow) b))))) X ((math-zerop X (setq q (math-sub (math-mul 4 (math-mul a c)) (math-sqr b)))) X (let ((part (math-div b (math-mul 2 c)))) X (math-mul-thru (math-pow c vpow) X (math-integral-q02 part 1 nil X (math-add math-integ-var part) X (* vpow 2))))) X ((progn X (setq part (math-add (math-mul 2 (math-mul c math-integ-var)) b)) X (> vpow 1)) X (let ((n (1- vpow))) X (math-add (math-div part (math-mul n (math-mul q (math-pow v n)))) X (math-mul-thru (math-div (math-mul (- (* 4 n) 2) c) X (math-mul n q)) X (math-integral-q02 a b c v n))))) X ((math-guess-if-neg q) X (setq rq (list 'calcFunc-sqrt (math-neg q))) X ;;(math-div-thru (list 'calcFunc-ln X ;; (math-div (math-sub part rq) X ;; (math-add part rq))) X ;; rq) X (math-div (math-mul -2 (list 'calcFunc-arctanh X (math-div part rq))) X rq)) X (t X (setq rq (list 'calcFunc-sqrt q)) X (math-div (math-mul 2 (math-to-radians-2 X (list 'calcFunc-arctan X (math-div part rq)))) X rq)))) ) X X X X (defun calcFunc-table (expr var &optional low high step) X (or low (setq low '(neg (var inf var-inf)) high '(var inf var-inf))) X (or high (setq high low low 1)) X (and (or (math-infinitep low) (math-infinitep high)) X (not step) X (math-scan-for-limits expr)) X (and step (math-zerop step) (math-reject-arg step 'nonzerop)) X (let ((known (+ (if (Math-objectp low) 1 0) X (if (Math-objectp high) 1 0) X (if (or (null step) (Math-objectp step)) 1 0))) X (count '(var inf var-inf)) X vec) X (or (= known 2) ; handy optimization X (equal high '(var inf var-inf)) X (progn X (setq count (math-div (math-sub high low) (or step 1))) X (or (Math-objectp count) X (setq count (math-simplify count))) X (if (Math-messy-integerp count) X (setq count (math-trunc count))))) X (if (Math-negp count) X (setq count -1)) X (if (integerp count) X (let ((var-DUMMY nil) X (vec math-tabulate-initial) X (math-working-step-2 (1+ count)) X (math-working-step 0)) X (setq expr (math-evaluate-expr X (math-expr-subst expr var '(var DUMMY var-DUMMY)))) X (while (>= count 0) X (setq math-working-step (1+ math-working-step) X var-DUMMY low X vec (cond ((eq math-tabulate-function 'calcFunc-sum) X (math-add vec (math-evaluate-expr expr))) X ((eq math-tabulate-function 'calcFunc-prod) X (math-mul vec (math-evaluate-expr expr))) X (t X (cons (math-evaluate-expr expr) vec))) X low (math-add low (or step 1)) X count (1- count))) X (if math-tabulate-function X vec X (cons 'vec (nreverse vec)))) X (if (Math-integerp count) X (calc-record-why 'fixnump high) X (if (Math-num-integerp low) X (if (Math-num-integerp high) X (calc-record-why 'integerp step) X (calc-record-why 'integerp high)) X (calc-record-why 'integerp low))) X (append (list (or math-tabulate-function 'calcFunc-table) X expr var) X (and (not (and (equal low '(neg (var inf var-inf))) X (equal high '(var inf var-inf)))) X (list low high)) X (and step (list step))))) ) X (setq math-tabulate-initial nil) (setq math-tabulate-function nil) X (defun math-scan-for-limits (x) X (cond ((Math-primp x)) X ((and (eq (car x) 'calcFunc-subscr) X (Math-vectorp (nth 1 x)) X (math-expr-contains (nth 2 x) var)) X (let* ((calc-next-why nil) X (low-val (math-solve-for (nth 2 x) 1 var nil)) X (high-val (math-solve-for (nth 2 x) (1- (length (nth 1 x))) X var nil)) X temp) X (and low-val (math-realp low-val) X high-val (math-realp high-val)) X (and (Math-lessp high-val low-val) X (setq temp low-val low-val high-val high-val temp)) X (setq low (math-max low (math-ceiling low-val)) X high (math-min high (math-floor high-val))))) X (t X (while (setq x (cdr x)) X (math-scan-for-limits (car x))))) ) X X (defun calcFunc-sum (expr var &optional low high step) X (if math-disable-sums (math-reject-arg)) X (let* ((res (let* ((calc-internal-prec (+ calc-internal-prec 2))) X (math-sum-rec expr var low high step))) X (math-disable-sums t)) X (math-normalize res)) ) (setq math-disable-sums nil) X (defun math-sum-rec (expr var &optional low high step) X (or low (setq low '(neg (var inf var-inf)) high '(var inf var-inf))) X (and low (not high) (setq high low low 1)) X (let (t1 t2 val) X (setq val X (cond X ((not (math-expr-contains expr var)) X (math-mul expr (math-add (math-div (math-sub high low) (or step 1)) X 1))) X ((and step (not (math-equal-int step 1))) X (if (math-negp step) X (math-sum-rec expr var high low (math-neg step)) X (let ((lo (math-simplify (math-div low step)))) X (if (math-known-num-integerp lo) X (math-sum-rec (math-normalize X (math-expr-subst expr var X (math-mul step var))) X var lo (math-simplify (math-div high step))) X (math-sum-rec (math-normalize X (math-expr-subst expr var X (math-add (math-mul step var) X low))) X var 0 X (math-simplify (math-div (math-sub high low) X step))))))) X ((memq (setq t1 (math-compare low high)) '(0 1)) X (if (eq t1 0) X (math-expr-subst expr var low) X 0)) X ((setq t1 (math-is-polynomial expr var 20)) X (let ((poly nil) X (n 0)) X (while t1 X (setq poly (math-poly-mix poly 1 X (math-sum-integer-power n) (car t1)) X n (1+ n) X t1 (cdr t1))) X (setq n (math-build-polynomial-expr poly high)) X (if (memq low '(0 1)) X n X (math-sub n (math-build-polynomial-expr poly X (math-sub low 1)))))) X ((and (memq (car expr) '(+ -)) X (setq t1 (math-sum-rec (nth 1 expr) var low high) X t2 (math-sum-rec (nth 2 expr) var low high)) X (not (and (math-expr-calls t1 '(calcFunc-sum)) X (math-expr-calls t2 '(calcFunc-sum))))) X (list (car expr) t1 t2)) X ((and (eq (car expr) '*) X (setq t1 (math-sum-const-factors expr var))) X (math-mul (car t1) (math-sum-rec (cdr t1) var low high))) X ((and (eq (car expr) '*) (memq (car-safe (nth 1 expr)) '(+ -))) X (math-sum-rec (math-add-or-sub (math-mul (nth 1 (nth 1 expr)) X (nth 2 expr)) X (math-mul (nth 2 (nth 1 expr)) X (nth 2 expr)) X nil (eq (car (nth 1 expr)) '-)) X var low high)) X ((and (eq (car expr) '*) (memq (car-safe (nth 2 expr)) '(+ -))) X (math-sum-rec (math-add-or-sub (math-mul (nth 1 expr) X (nth 1 (nth 2 expr))) X (math-mul (nth 1 expr) X (nth 2 (nth 2 expr))) X nil (eq (car (nth 2 expr)) '-)) X var low high)) X ((and (eq (car expr) '/) X (not (math-primp (nth 1 expr))) X (setq t1 (math-sum-const-factors (nth 1 expr) var))) X (math-mul (car t1) X (math-sum-rec (math-div (cdr t1) (nth 2 expr)) X var low high))) X ((and (eq (car expr) '/) X (setq t1 (math-sum-const-factors (nth 2 expr) var))) X (math-div (math-sum-rec (math-div (nth 1 expr) (cdr t1)) X var low high) X (car t1))) X ((eq (car expr) 'neg) X (math-neg (math-sum-rec (nth 1 expr) var low high))) X ((and (eq (car expr) '^) X (not (math-expr-contains (nth 1 expr) var)) X (setq t1 (math-is-polynomial (nth 2 expr) var 1))) X (let ((x (math-pow (nth 1 expr) (nth 1 t1)))) X (math-div (math-mul (math-sub (math-pow x (math-add 1 high)) X (math-pow x low)) X (math-pow (nth 1 expr) (car t1))) X (math-sub x 1)))) X ((and (setq t1 (math-to-exponentials expr)) X (setq t1 (math-sum-rec t1 var low high)) X (not (math-expr-calls t1 '(calcFunc-sum)))) X (math-to-exps t1)) X ((memq (car expr) '(calcFunc-ln calcFunc-log10)) X (list (car expr) (calcFunc-prod (nth 1 expr) var low high))) X ((and (eq (car expr) 'calcFunc-log) X (= (length expr) 3) X (not (math-expr-contains (nth 2 expr) var))) X (list 'calcFunc-log X (calcFunc-prod (nth 1 expr) var low high) X (nth 2 expr))))) X (if (equal val '(var nan var-nan)) (setq val nil)) X (or val X (let* ((math-tabulate-initial 0) X (math-tabulate-function 'calcFunc-sum)) X (calcFunc-table expr var low high)))) ) X (defun calcFunc-asum (expr var low &optional high step no-mul-flag) X (or high (setq high low low 1)) X (if (and step (not (math-equal-int step 1))) X (if (math-negp step) X (math-mul (math-pow -1 low) X (calcFunc-asum expr var high low (math-neg step) t)) X (let ((lo (math-simplify (math-div low step)))) X (if (math-num-integerp lo) X (calcFunc-asum (math-normalize X (math-expr-subst expr var X (math-mul step var))) X var lo (math-simplify (math-div high step))) X (calcFunc-asum (math-normalize X (math-expr-subst expr var X (math-add (math-mul step var) X low))) X var 0 X (math-simplify (math-div (math-sub high low) X step)))))) X (math-mul (if no-mul-flag 1 (math-pow -1 low)) X (calcFunc-sum (math-mul (math-pow -1 var) expr) var low high))) ) X (defun math-sum-const-factors (expr var) X (let ((const nil) X (not-const nil) X (p expr)) X (while (eq (car-safe p) '*) X (if (math-expr-contains (nth 1 p) var) X (setq not-const (cons (nth 1 p) not-const)) X (setq const (cons (nth 1 p) const))) X (setq p (nth 2 p))) X (if (math-expr-contains p var) X (setq not-const (cons p not-const)) X (setq const (cons p const))) X (and const X (cons (let ((temp (car const))) X (while (setq const (cdr const)) X (setq temp (list '* (car const) temp))) X temp) X (let ((temp (or (car not-const) 1))) X (while (setq not-const (cdr not-const)) X (setq temp (list '* (car not-const) temp))) X temp)))) ) X ;; Following is from CRC Math Tables, 27th ed, pp. 52-53. (defun math-sum-integer-power (pow) X (let ((calc-prefer-frac t) X (n (length math-sum-int-pow-cache))) X (while (<= n pow) X (let* ((new (list 0 0)) X (lin new) X (pp (cdr (nth (1- n) math-sum-int-pow-cache))) X (p 2) X (sum 0) X q) X (while pp X (setq q (math-div (car pp) p) X new (cons (math-mul q n) new) X sum (math-add sum q) X p (1+ p) X pp (cdr pp))) X (setcar lin (math-sub 1 (math-mul n sum))) X (setq math-sum-int-pow-cache X (nconc math-sum-int-pow-cache (list (nreverse new))) X n (1+ n)))) X (nth pow math-sum-int-pow-cache)) ) (setq math-sum-int-pow-cache (list '(0 1))) X (defun math-to-exponentials (expr) X (and (consp expr) X (= (length expr) 2) X (let ((x (nth 1 expr)) X (pi (if calc-symbolic-mode '(var pi var-pi) (math-pi))) X (i (if calc-symbolic-mode '(var i var-i) '(cplx 0 1)))) X (cond ((eq (car expr) 'calcFunc-exp) X (list '^ '(var e var-e) x)) X ((eq (car expr) 'calcFunc-sin) X (or (eq calc-angle-mode 'rad) X (setq x (list '/ (list '* x pi) 180))) X (list '/ (list '- X (list '^ '(var e var-e) (list '* x i)) X (list '^ '(var e var-e) X (list 'neg (list '* x i)))) X (list '* 2 i))) X ((eq (car expr) 'calcFunc-cos) X (or (eq calc-angle-mode 'rad) X (setq x (list '/ (list '* x pi) 180))) X (list '/ (list '+ X (list '^ '(var e var-e) X (list '* x i)) X (list '^ '(var e var-e) X (list 'neg (list '* x i)))) X 2)) X ((eq (car expr) 'calcFunc-sinh) X (list '/ (list '- X (list '^ '(var e var-e) x) X (list '^ '(var e var-e) (list 'neg x))) X 2)) X ((eq (car expr) 'calcFunc-cosh) X (list '/ (list '+ X (list '^ '(var e var-e) x) X (list '^ '(var e var-e) (list 'neg x))) X 2)) X (t nil)))) ) X (defun math-to-exps (expr) X (cond (calc-symbolic-mode expr) X ((Math-primp expr) X (if (equal expr '(var e var-e)) (math-e) expr)) X ((and (eq (car expr) '^) X (equal (nth 1 expr) '(var e var-e))) X (list 'calcFunc-exp (nth 2 expr))) X (t X (cons (car expr) (mapcar 'math-to-exps (cdr expr))))) ) X X (defun calcFunc-prod (expr var &optional low high step) X (if math-disable-prods (math-reject-arg)) X (let* ((res (let* ((calc-internal-prec (+ calc-internal-prec 2))) X (math-prod-rec expr var low high step))) X (math-disable-prods t)) X (math-normalize res)) ) (setq math-disable-prods nil) X (defun math-prod-rec (expr var &optional low high step) X (or low (setq low '(neg (var inf var-inf)) high '(var inf var-inf))) X (and low (not high) (setq high '(var inf var-inf))) X (let (t1 t2 t3 val) X (setq val X (cond X ((not (math-expr-contains expr var)) X (math-pow expr (math-add (math-div (math-sub high low) (or step 1)) X 1))) X ((and step (not (math-equal-int step 1))) X (if (math-negp step) X (math-prod-rec expr var high low (math-neg step)) X (let ((lo (math-simplify (math-div low step)))) X (if (math-known-num-integerp lo) X (math-prod-rec (math-normalize X (math-expr-subst expr var X (math-mul step var))) X var lo (math-simplify (math-div high step))) X (math-prod-rec (math-normalize X (math-expr-subst expr var X (math-add (math-mul step X var) X low))) X var 0 X (math-simplify (math-div (math-sub high low) X step))))))) X ((and (memq (car expr) '(* /)) X (setq t1 (math-prod-rec (nth 1 expr) var low high) X t2 (math-prod-rec (nth 2 expr) var low high)) X (not (and (math-expr-calls t1 '(calcFunc-prod)) X (math-expr-calls t2 '(calcFunc-prod))))) X (list (car expr) t1 t2)) X ((and (eq (car expr) '^) X (not (math-expr-contains (nth 2 expr) var))) X (math-pow (math-prod-rec (nth 1 expr) var low high) X (nth 2 expr))) X ((and (eq (car expr) '^) X (not (math-expr-contains (nth 1 expr) var))) X (math-pow (nth 1 expr) X (calcFunc-sum (nth 2 expr) var low high))) X ((eq (car expr) 'sqrt) X (math-normalize (list 'calcFunc-sqrt X (list 'calcFunc-prod (nth 1 expr) X var low high)))) X ((eq (car expr) 'neg) X (math-mul (math-pow -1 (math-add (math-sub high low) 1)) X (math-prod-rec (nth 1 expr) var low high))) X ((eq (car expr) 'calcFunc-exp) X (list 'calcFunc-exp (calcFunc-sum (nth 1 expr) var low high))) X ((and (setq t1 (math-is-polynomial expr var 1)) X (setq t2 X (cond X ((or (and (math-equal-int (nth 1 t1) 1) X (setq low (math-simplify X (math-add low (car t1))) X high (math-simplify X (math-add high (car t1))))) X (and (math-equal-int (nth 1 t1) -1) X (setq t2 low X low (math-simplify X (math-sub (car t1) high)) X high (math-simplify X (math-sub (car t1) t2))))) X (if (or (math-zerop low) (math-zerop high)) X 0 X (if (and (or (math-negp low) (math-negp high)) X (or (math-num-integerp low) X (math-num-integerp high))) X (if (math-posp high) X 0 X (math-mul (math-pow -1 X (math-add X (math-add low high) 1)) X (list '/ X (list 'calcFunc-fact X (math-neg low)) X (list 'calcFunc-fact X (math-sub -1 high))))) X (list '/ X (list 'calcFunc-fact high) X (list 'calcFunc-fact (math-sub low 1)))))) X ((and (or (and (math-equal-int (nth 1 t1) 2) X (setq t2 (math-simplify X (math-add (math-mul low 2) X (car t1))) X t3 (math-simplify X (math-add (math-mul high 2) X (car t1))))) X (and (math-equal-int (nth 1 t1) -2) X (setq t2 (math-simplify X (math-sub (car t1) X (math-mul high 2))) X t3 (math-simplify X (math-sub (car t1) X (math-mul low X 2)))))) X (or (math-integerp t2) X (and (math-messy-integerp t2) X (setq t2 (math-trunc t2))) X (math-integerp t3) X (and (math-messy-integerp t3) X (setq t3 (math-trunc t3))))) X (if (or (math-zerop t2) (math-zerop t3)) X 0 X (if (or (math-evenp t2) (math-evenp t3)) X (if (or (math-negp t2) (math-negp t3)) X (if (math-posp high) X 0 X (list '/ X (list 'calcFunc-dfact X (math-neg t2)) X (list 'calcFunc-dfact X (math-sub -2 t3)))) X (list '/ X (list 'calcFunc-dfact t3) X (list 'calcFunc-dfact X (math-sub t2 2)))) X (if (math-negp t3) X (list '* X (list '^ -1 X (list '/ (list '- (list '- t2 t3) X 2) X 2)) X (list '/ X (list 'calcFunc-dfact X (math-neg t2)) X (list 'calcFunc-dfact X (math-sub -2 t3)))) X (if (math-posp t2) X (list '/ X (list 'calcFunc-dfact t3) X (list 'calcFunc-dfact X (math-sub t2 2))) X nil)))))))) X t2))) X (if (equal val '(var nan var-nan)) (setq val nil)) X (or val X (let* ((math-tabulate-initial 1) X (math-tabulate-function 'calcFunc-prod)) X (calcFunc-table expr var low high)))) ) X X X X ;;; Attempt to reduce lhs = rhs to solve-var = rhs', where solve-var appears ;;; in lhs but not in rhs or rhs'; return rhs'. ;;; Uses global values: solve-*. (defun math-try-solve-for (lhs rhs &optional sign no-poly) X (let (t1 t2 t3) X (cond ((equal lhs solve-var) X (setq math-solve-sign sign) X (if (eq solve-full 'all) X (let ((vec (list 'vec (math-evaluate-expr rhs))) X newvec var p) X (while math-solve-ranges X (setq p (car math-solve-ranges) X var (car p) X newvec (list 'vec)) X (while (setq p (cdr p)) X (setq newvec (nconc newvec X (cdr (math-expr-subst X vec var (car p)))))) X (setq vec newvec X math-solve-ranges (cdr math-solve-ranges))) X (math-normalize vec)) X rhs)) X ((Math-primp lhs) X nil) X ((and (eq (car lhs) '-) X (eq (car-safe (nth 1 lhs)) (car-safe (nth 2 lhs))) X (Math-zerop rhs) X (= (length (nth 1 lhs)) 2) X (= (length (nth 2 lhs)) 2) X (setq t1 (get (car (nth 1 lhs)) 'math-inverse)) X (setq t2 (funcall t1 '(var SOLVEDUM SOLVEDUM))) X (eq (math-expr-contains-count t2 '(var SOLVEDUM SOLVEDUM)) 1) X (setq t3 (math-solve-above-dummy t2)) X (setq t1 (math-try-solve-for (math-sub (nth 1 (nth 1 lhs)) X (math-expr-subst X t2 t3 X (nth 1 (nth 2 lhs)))) X 0))) X t1) X ((eq (car lhs) 'neg) X (math-try-solve-for (nth 1 lhs) (math-neg rhs) X (and sign (- sign)))) X ((and (not (eq solve-full 't)) (math-try-solve-prod))) X ((and (not no-poly) X (setq t2 (math-decompose-poly lhs solve-var 15 rhs))) X (setq t1 (cdr (nth 1 t2)) X t1 (let ((math-solve-ranges math-solve-ranges)) X (cond ((= (length t1) 5) X (apply 'math-solve-quartic (car t2) t1)) X ((= (length t1) 4) X (apply 'math-solve-cubic (car t2) t1)) X ((= (length t1) 3) X (apply 'math-solve-quadratic (car t2) t1)) X ((= (length t1) 2) X (apply 'math-solve-linear (car t2) sign t1)) X (solve-full X (math-poly-all-roots (car t2) t1)) X (calc-symbolic-mode nil) X (t X (math-try-solve-for X (car t2) X (math-poly-any-root (reverse t1) 0 t) X nil t))))) X (if t1 X (if (eq (nth 2 t2) 1) X t1 X (math-solve-prod t1 (math-try-solve-for (nth 2 t2) 0 nil t))) X (calc-record-why "*Unable to find a symbolic solution") X nil)) X ((and (math-solve-find-root-term lhs nil) X (eq (math-expr-contains-count lhs t1) 1)) ; just in case X (math-try-solve-for (math-simplify X (math-sub (if (or t3 (math-evenp t2)) X (math-pow t1 t2) X (math-neg (math-pow t1 t2))) X (math-expand-power X (math-sub (math-normalize X (math-expr-subst X lhs t1 0)) X rhs) X t2 solve-var))) X 0)) X ((eq (car lhs) '+) X (cond ((not (math-expr-contains (nth 1 lhs) solve-var)) X (math-try-solve-for (nth 2 lhs) X (math-sub rhs (nth 1 lhs)) X sign)) X ((not (math-expr-contains (nth 2 lhs) solve-var)) X (math-try-solve-for (nth 1 lhs) X (math-sub rhs (nth 2 lhs)) X sign)))) X ((eq (car lhs) 'calcFunc-eq) X (math-try-solve-for (math-sub (nth 1 lhs) (nth 2 lhs)) X rhs sign no-poly)) X ((eq (car lhs) '-) X (cond ((or (and (eq (car-safe (nth 1 lhs)) 'calcFunc-sin) X (eq (car-safe (nth 2 lhs)) 'calcFunc-cos)) X (and (eq (car-safe (nth 1 lhs)) 'calcFunc-cos) X (eq (car-safe (nth 2 lhs)) 'calcFunc-sin))) X (math-try-solve-for (math-sub (nth 1 lhs) X (list (car (nth 1 lhs)) X (math-sub X (math-quarter-circle t) X (nth 1 (nth 2 lhs))))) X rhs)) X ((not (math-expr-contains (nth 1 lhs) solve-var)) X (math-try-solve-for (nth 2 lhs) X (math-sub (nth 1 lhs) rhs) X (and sign (- sign)))) X ((not (math-expr-contains (nth 2 lhs) solve-var)) X (math-try-solve-for (nth 1 lhs) X (math-add rhs (nth 2 lhs)) X sign)))) X ((and (eq solve-full 't) (math-try-solve-prod))) X ((and (eq (car lhs) '%) X (not (math-expr-contains (nth 2 lhs) solve-var))) X (math-try-solve-for (nth 1 lhs) (math-add rhs X (math-solve-get-int X (nth 2 lhs))))) X ((eq (car lhs) 'calcFunc-log) X (cond ((not (math-expr-contains (nth 2 lhs) solve-var)) X (math-try-solve-for (nth 1 lhs) (math-pow (nth 2 lhs) rhs))) X ((not (math-expr-contains (nth 1 lhs) solve-var)) X (math-try-solve-for (nth 2 lhs) (math-pow X (nth 1 lhs) X (math-div 1 rhs)))))) X ((and (= (length lhs) 2) X (symbolp (car lhs)) X (setq t1 (get (car lhs) 'math-inverse)) X (setq t2 (funcall t1 rhs))) X (setq t1 (get (car lhs) 'math-inverse-sign)) X (math-try-solve-for (nth 1 lhs) (math-normalize t2) X (and sign t1 X (if (integerp t1) X (* t1 sign) X (funcall t1 lhs sign))))) X ((and (symbolp (car lhs)) X (setq t1 (get (car lhs) 'math-inverse-n)) X (setq t2 (funcall t1 lhs rhs))) X t2) X ((setq t1 (math-expand-formula lhs)) X (math-try-solve-for t1 rhs sign)) X (t X (calc-record-why "*No inverse known" lhs) X nil))) ) X (setq math-solve-ranges nil) X (defun math-try-solve-prod () X (cond ((eq (car lhs) '*) X (cond ((not (math-expr-contains (nth 1 lhs) solve-var)) X (math-try-solve-for (nth 2 lhs) X (math-div rhs (nth 1 lhs)) X (math-solve-sign sign (nth 1 lhs)))) X ((not (math-expr-contains (nth 2 lhs) solve-var)) X (math-try-solve-for (nth 1 lhs) X (math-div rhs (nth 2 lhs)) X (math-solve-sign sign (nth 2 lhs)))) X ((Math-zerop rhs) X (math-solve-prod (let ((math-solve-ranges math-solve-ranges)) X (math-try-solve-for (nth 2 lhs) 0)) X (math-try-solve-for (nth 1 lhs) 0))))) X ((eq (car lhs) '/) X (cond ((not (math-expr-contains (nth 1 lhs) solve-var)) X (math-try-solve-for (nth 2 lhs) X (math-div (nth 1 lhs) rhs) X (math-solve-sign sign (nth 1 lhs)))) X ((not (math-expr-contains (nth 2 lhs) solve-var)) X (math-try-solve-for (nth 1 lhs) X (math-mul rhs (nth 2 lhs)) X (math-solve-sign sign (nth 2 lhs)))) X ((setq t1 (math-try-solve-for (math-sub (nth 1 lhs) X (math-mul (nth 2 lhs) X rhs)) X 0)) X t1))) X ((eq (car lhs) '^) X (cond ((not (math-expr-contains (nth 1 lhs) solve-var)) X (math-try-solve-for X (nth 2 lhs) X (math-add (math-normalize X (list 'calcFunc-log rhs (nth 1 lhs))) X (math-div X (math-mul 2 X (math-mul '(var pi var-pi) X (math-solve-get-int X '(var i var-i)))) X (math-normalize X (list 'calcFunc-ln (nth 1 lhs))))))) X ((not (math-expr-contains (nth 2 lhs) solve-var)) X (cond ((and (integerp (nth 2 lhs)) X (>= (nth 2 lhs) 2) X (setq t1 (math-integer-log2 (nth 2 lhs)))) X (setq t2 rhs) X (if (and (eq solve-full t) X (math-known-realp (nth 1 lhs))) X (progn X (while (>= (setq t1 (1- t1)) 0) X (setq t2 (list 'calcFunc-sqrt t2))) X (setq t2 (math-solve-get-sign t2))) X (while (>= (setq t1 (1- t1)) 0) X (setq t2 (math-solve-get-sign X (math-normalize X (list 'calcFunc-sqrt t2)))))) X (math-try-solve-for X (nth 1 lhs) X (math-normalize t2))) X ((math-looks-negp (nth 2 lhs)) X (math-try-solve-for X (list '^ (nth 1 lhs) (math-neg (nth 2 lhs))) X (math-div 1 rhs))) X ((and (eq solve-full t) X (Math-integerp (nth 2 lhs)) X (math-known-realp (nth 1 lhs))) X (setq t1 (math-normalize X (list 'calcFunc-nroot rhs (nth 2 lhs)))) X (if (math-evenp (nth 2 lhs)) X (setq t1 (math-solve-get-sign t1))) X (math-try-solve-for X (nth 1 lhs) t1 X (and sign X (math-oddp (nth 2 lhs)) X (math-solve-sign sign (nth 2 lhs))))) X (t (math-try-solve-for X (nth 1 lhs) X (math-mul X (math-normalize X (list 'calcFunc-exp X (if (Math-realp (nth 2 lhs)) X (math-div (math-mul X '(var pi var-pi) X (math-solve-get-int X '(var i var-i) X (and (integerp (nth 2 lhs)) X (math-abs X (nth 2 lhs))))) X (math-div (nth 2 lhs) 2)) X (math-div (math-mul X 2 X (math-mul X '(var pi var-pi) X (math-solve-get-int X '(var i var-i) X (and (integerp (nth 2 lhs)) X (math-abs X (nth 2 lhs)))))) X (nth 2 lhs))))) X (math-normalize X (list 'calcFunc-nroot X rhs X (nth 2 lhs)))) X (and sign X (math-oddp (nth 2 lhs)) X (math-solve-sign sign (nth 2 lhs))))))))) X (t nil)) ) X (defun math-solve-prod (lsoln rsoln) X (cond ((null lsoln) X rsoln) X ((null rsoln) X lsoln) X ((eq solve-full 'all) X (cons 'vec (append (cdr lsoln) (cdr rsoln)))) X (solve-full X (list 'calcFunc-if X (list 'calcFunc-gt (math-solve-get-sign 1) 0) X lsoln X rsoln)) X (t lsoln)) ) X ;;; This deals with negative, fractional, and symbolic powers of "x". (defun math-solve-poly-funny-powers (sub-rhs) ; uses "t1", "t2" X (setq t1 lhs) X (let ((pp math-poly-neg-powers) X fac) X (while pp X (setq fac (math-pow (car pp) (or math-poly-mult-powers 1)) X t1 (math-mul t1 fac) X rhs (math-mul rhs fac) X pp (cdr pp)))) X (if sub-rhs (setq t1 (math-sub t1 rhs))) X (let ((math-poly-neg-powers nil)) X (setq t2 (math-mul (or math-poly-mult-powers 1) X (let ((calc-prefer-frac t)) X (math-div 1 math-poly-frac-powers))) X t1 (math-is-polynomial (math-simplify (calcFunc-expand t1)) b 50))) ) X ;;; This converts "a x^8 + b x^5 + c x^2" to "(a (x^3)^2 + b (x^3) + c) * x^2". (defun math-solve-crunch-poly (max-degree) ; uses "t1", "t3" X (let ((count 0)) X (while (and t1 (Math-zerop (car t1))) X (setq t1 (cdr t1) X count (1+ count))) X (and t1 X (let* ((degree (1- (length t1))) X (scale degree)) X (while (and (> scale 1) (= (car t3) 1)) X (and (= (% degree scale) 0) X (let ((p t1) X (n 0) X (new-t1 nil) X (okay t)) X (while (and p okay) X (if (= (% n scale) 0) X (setq new-t1 (nconc new-t1 (list (car p)))) X (or (Math-zerop (car p)) X (setq okay nil))) X (setq p (cdr p) X n (1+ n))) X (if okay X (setq t3 (cons scale (cdr t3)) X t1 new-t1)))) X (setq scale (1- scale))) X (setq t3 (list (math-mul (car t3) t2) (math-mul count t2))) X (<= (1- (length t1)) max-degree)))) ) X (defun calcFunc-poly (expr var &optional degree) X (if degree X (or (natnump degree) (math-reject-arg degree 'fixnatnump)) X (setq degree 50)) X (let ((p (math-is-polynomial expr var degree 'gen))) X (if p X (if (equal p '(0)) X (list 'vec) X (cons 'vec p)) X (math-reject-arg expr "Expected a polynomial"))) ) X (defun calcFunc-gpoly (expr var &optional degree) X (if degree X (or (natnump degree) (math-reject-arg degree 'fixnatnump)) X (setq degree 50)) X (let* ((math-poly-base-variable var) X (d (math-decompose-poly expr var degree nil))) X (if d X (cons 'vec d) X (math-reject-arg expr "Expected a polynomial"))) ) X (defun math-decompose-poly (lhs solve-var degree sub-rhs) X (let ((rhs (or sub-rhs 1)) X t1 t2 t3) X (setq t2 (math-polynomial-base X lhs X (function X (lambda (b) X (let ((math-poly-neg-powers '(1)) X (math-poly-mult-powers nil) X (math-poly-frac-powers 1) X (math-poly-exp-base t)) X (and (not (equal b lhs)) X (or (not (memq (car-safe b) '(+ -))) sub-rhs) X (setq t3 '(1 0) t2 1 X t1 (math-is-polynomial lhs b 50)) X (if (and (equal math-poly-neg-powers '(1)) X (memq math-poly-mult-powers '(nil 1)) X (eq math-poly-frac-powers 1) X sub-rhs) X (setq t1 (cons (math-sub (car t1) rhs) X (cdr t1))) X (math-solve-poly-funny-powers sub-rhs)) X (math-solve-crunch-poly degree) X (or (math-expr-contains b solve-var) X (math-expr-contains (car t3) solve-var)))))))) X (if t2 X (list (math-pow t2 (car t3)) X (cons 'vec t1) X (if sub-rhs X (math-pow t2 (nth 1 t3)) X (math-div (math-pow t2 (nth 1 t3)) rhs))))) ) X (defun math-solve-linear (var sign b a) X (math-try-solve-for var X (math-div (math-neg b) a) X (math-solve-sign sign a) X t) ) X (defun math-solve-quadratic (var c b a) X (math-try-solve-for X var X (if (math-looks-evenp b) X (let ((halfb (math-div b 2))) X (math-div X (math-add X (math-neg halfb) X (math-solve-get-sign X (math-normalize X (list 'calcFunc-sqrt X (math-add (math-sqr halfb) X (math-mul (math-neg c) a)))))) X a)) X (math-div X (math-add X (math-neg b) X (math-solve-get-sign X (math-normalize X (list 'calcFunc-sqrt X (math-add (math-sqr b) X (math-mul 4 (math-mul (math-neg c) a))))))) X (math-mul 2 a))) X nil t) ) X (defun math-solve-cubic (var d c b a) X (let* ((p (math-div b a)) X (q (math-div c a)) X (r (math-div d a)) X (psqr (math-sqr p)) X (aa (math-sub q (math-div psqr 3))) X (bb (math-add r X (math-div (math-sub (math-mul 2 (math-mul psqr p)) X (math-mul 9 (math-mul p q))) X 27))) X m) X (if (Math-zerop aa) X (math-try-solve-for (math-pow (math-add var (math-div p 3)) 3) X (math-neg bb) nil t) X (if (Math-zerop bb) X (math-try-solve-for X (math-mul (math-add var (math-div p 3)) X (math-add (math-sqr (math-add var (math-div p 3))) X aa)) X 0 nil t) X (setq m (math-mul 2 (list 'calcFunc-sqrt (math-div aa -3)))) X (math-try-solve-for X var X (math-sub X (math-normalize X (math-mul X m X (list 'calcFunc-cos X (math-div X (math-sub (list 'calcFunc-arccos X (math-div (math-mul 3 bb) X (math-mul aa m))) X (math-mul 2 X (math-mul X (math-add 1 (math-solve-get-int X 1 3)) X (math-half-circle X calc-symbolic-mode)))) X 3)))) X (math-div p 3)) X nil t)))) ) X (defun math-solve-quartic (var d c b a aa) X (setq a (math-div a aa)) X (setq b (math-div b aa)) X (setq c (math-div c aa)) X (setq d (math-div d aa)) X (math-try-solve-for X var X (let* ((asqr (math-sqr a)) X (asqr4 (math-div asqr 4)) X (y (let ((solve-full nil) X calc-next-why) X (math-solve-cubic solve-var X (math-sub (math-sub X (math-mul 4 (math-mul b d)) X (math-mul asqr d)) X (math-sqr c)) X (math-sub (math-mul a c) X (math-mul 4 d)) X (math-neg b) X 1))) X (rsqr (math-add (math-sub asqr4 b) y)) X (r (list 'calcFunc-sqrt rsqr)) X (sign1 (math-solve-get-sign 1)) X (de (list 'calcFunc-sqrt X (math-add X (math-sub (math-mul 3 asqr4) X (math-mul 2 b)) X (if (Math-zerop rsqr) X (math-mul X 2 X (math-mul sign1 X (list 'calcFunc-sqrt X (math-sub (math-sqr y) X (math-mul 4 d))))) X (math-sub X (math-mul sign1 X (math-div X (math-sub (math-sub X (math-mul 4 (math-mul a b)) X (math-mul 8 c)) X (math-mul asqr a)) X (math-mul 4 r))) X rsqr)))))) X (math-normalize X (math-sub (math-add (math-mul sign1 (math-div r 2)) X (math-solve-get-sign (math-div de 2))) X (math-div a 4)))) X nil t) ) X (defun math-poly-all-roots (var p &optional math-factoring) X (catch 'ouch X (let* ((math-symbolic-solve calc-symbolic-mode) X (roots nil) X (deg (1- (length p))) X (orig-p (reverse p)) X (math-int-coefs nil) X (math-int-scale nil) X (math-double-roots nil) X (math-int-factors nil) X (math-int-threshold nil) X (pp p)) X ;; If rational coefficients, look for exact rational factors. X (while (and pp (Math-ratp (car pp))) X (setq pp (cdr pp))) X (if pp X (if (or math-factoring math-symbolic-solve) X (throw 'ouch nil)) X (let ((lead (car orig-p)) X (calc-prefer-frac t) X (scale (apply 'math-lcm-denoms p))) X (setq math-int-scale (math-abs (math-mul scale lead)) X math-int-threshold (math-div '(float 5 -2) math-int-scale) X math-int-coefs (cdr (math-div (cons 'vec orig-p) lead))))) X (if (> deg 4) X (let ((calc-prefer-frac nil) X (calc-symbolic-mode nil) X (pp p) X (def-p (copy-sequence orig-p))) X (while pp X (if (Math-numberp (car pp)) X (setq pp (cdr pp)) X (throw 'ouch nil))) X (while (> deg (if math-symbolic-solve 2 4)) X (let* ((x (math-poly-any-root def-p '(float 0 0) nil)) X b c pp) X (if (and (eq (car-safe x) 'cplx) X (math-nearly-zerop (nth 2 x) (nth 1 x))) X (setq x (calcFunc-re x))) X (or math-factoring X (setq roots (cons x roots))) X (or (math-numberp x) X (setq x (math-evaluate-expr x))) X (setq pp def-p X b (car def-p)) X (while (setq pp (cdr pp)) X (setq c (car pp)) X (setcar pp b) X (setq b (math-add (math-mul x b) c))) X (setq def-p (cdr def-p) X deg (1- deg)))) X (setq p (reverse def-p)))) X (if (> deg 1) X (let ((solve-var '(var DUMMY var-DUMMY)) X (math-solve-sign nil) X (math-solve-ranges nil) X (solve-full 'all)) X (if (= (length p) (length math-int-coefs)) X (setq p (reverse math-int-coefs))) X (setq roots (append (cdr (apply (cond ((= deg 2) X 'math-solve-quadratic) X ((= deg 3) X 'math-solve-cubic) X (t X 'math-solve-quartic)) X solve-var p)) X roots))) X (if (> deg 0) X (setq roots (cons (math-div (math-neg (car p)) (nth 1 p)) X roots)))) X (if math-factoring X (progn X (while roots X (math-poly-integer-root (car roots)) X (setq roots (cdr roots))) X (list math-int-factors (nreverse math-int-coefs) math-int-scale)) X (let ((vec nil) res) X (while roots X (let ((root (car roots)) X (solve-full (and solve-full 'all))) X (if (math-floatp root) X (setq root (math-poly-any-root orig-p root t))) X (setq vec (append vec X (cdr (or (math-try-solve-for var root nil t) X (throw 'ouch nil)))))) X (setq roots (cdr roots))) X (setq vec (cons 'vec (nreverse vec))) X (if math-symbolic-solve X (setq vec (math-normalize vec))) X (if (eq solve-full t) X (list 'calcFunc-subscr X vec X (math-solve-get-int 1 (1- (length orig-p)) 1)) X vec))))) ) (setq math-symbolic-solve nil) X (defun math-lcm-denoms (&rest fracs) X (let ((den 1)) X (while fracs X (if (eq (car-safe (car fracs)) 'frac) X (setq den (calcFunc-lcm den (nth 2 (car fracs))))) X (setq fracs (cdr fracs))) X den) ) X (defun math-poly-any-root (p x polish) ; p is a reverse poly coeff list X (let* ((newt (if (math-zerop x) X (math-poly-newton-root X p '(cplx (float 123 -6) (float 1 -4)) 4) X (math-poly-newton-root p x 4))) X (res (if (math-zerop (cdr newt)) X (car newt) X (if (and (math-lessp (cdr newt) '(float 1 -3)) (not polish)) X (setq newt (math-poly-newton-root p (car newt) 30))) X (if (math-zerop (cdr newt)) X (car newt) X (math-poly-laguerre-root p x polish))))) X (and math-symbolic-solve (math-floatp res) X (throw 'ouch nil)) X res) ) X (defun math-poly-newton-root (p x iters) X (let* ((calc-prefer-frac nil) X (calc-symbolic-mode nil) X (try-integer math-int-coefs) X (dx x) b d) X (while (and (> (setq iters (1- iters)) 0) X (let ((pp p)) X (math-working "newton" x) X (setq b (car p) X d 0) X (while (setq pp (cdr pp)) X (setq d (math-add (math-mul x d) b) X b (math-add (math-mul x b) (car pp)))) X (not (math-zerop d))) X (progn X (setq dx (math-div b d) X x (math-sub x dx)) X (if try-integer X (let ((adx (math-abs-approx dx))) X (and (math-lessp adx math-int-threshold) X (let ((iroot (math-poly-integer-root x))) X (if iroot X (setq x iroot dx 0) X (setq try-integer nil)))))) X (or (not (or (eq dx 0) X (math-nearly-zerop dx (math-abs-approx x)))) X (progn (setq dx 0) nil))))) X (cons x (if (math-zerop x) X 1 (math-div (math-abs-approx dx) (math-abs-approx x))))) ) X (defun math-poly-integer-root (x) X (and (math-lessp (calcFunc-xpon (math-abs-approx x)) calc-internal-prec) X math-int-coefs X (let* ((calc-prefer-frac t) X (xre (calcFunc-re x)) X (xim (calcFunc-im x)) X (xresq (math-sqr xre)) X (ximsq (math-sqr xim))) X (if (math-lessp ximsq (calcFunc-scf xresq -1)) X ;; Look for linear factor X (let* ((rnd (math-div (math-round (math-mul xre math-int-scale)) X math-int-scale)) X (icp math-int-coefs) X (rem (car icp)) X (newcoef nil)) X (while (setq icp (cdr icp)) X (setq newcoef (cons rem newcoef) X rem (math-add (car icp) X (math-mul rem rnd)))) X (and (math-zerop rem) X (progn X (setq math-int-coefs (nreverse newcoef) X math-int-factors (cons (list (math-neg rnd)) X math-int-factors)) X rnd))) X ;; Look for irreducible quadratic factor X (let* ((rnd1 (math-div (math-round X (math-mul xre (math-mul -2 math-int-scale))) X math-int-scale)) X (sqscale (math-sqr math-int-scale)) X (rnd0 (math-div (math-round (math-mul (math-add xresq ximsq) X sqscale)) X sqscale)) X (rem1 (car math-int-coefs)) X (icp (cdr math-int-coefs)) X (rem0 (car icp)) X (newcoef nil) X (found (assoc (list rnd0 rnd1 (math-posp xim)) X math-double-roots)) X this) X (if found X (setq math-double-roots (delq found math-double-roots) X rem0 0 rem1 0) X (while (setq icp (cdr icp)) X (setq this rem1 X newcoef (cons rem1 newcoef) X rem1 (math-sub rem0 (math-mul this rnd1)) X rem0 (math-sub (car icp) (math-mul this rnd0))))) X (and (math-zerop rem0) X (math-zerop rem1) X (let ((aa (math-div rnd1 -2))) X (or found (setq math-int-coefs (reverse newcoef) X math-double-roots (cons (list X (list X rnd0 rnd1 X (math-negp xim))) X math-double-roots) X math-int-factors (cons (cons rnd0 rnd1) X math-int-factors))) X (math-add aa X (let ((calc-symbolic-mode math-symbolic-solve)) X (math-mul (math-sqrt (math-sub (math-sqr aa) X rnd0)) X (if (math-negp xim) -1 1)))))))))) ) (setq math-int-coefs nil) X ;;; The following routine is from Numerical Recipes, section 9.5. (defun math-poly-laguerre-root (p x polish) X (let* ((calc-prefer-frac nil) X (calc-symbolic-mode nil) X (iters 0) X (m (1- (length p))) X (try-newt (not polish)) X (tried-newt nil) X b d f x1 dx dxold) X (while X (and (or (< (setq iters (1+ iters)) 50) X (math-reject-arg x "*Laguerre's method failed to converge")) X (let ((err (math-abs-approx (car p))) X (abx (math-abs-approx x)) X (pp p)) X (setq b (car p) X d 0 f 0) X (while (setq pp (cdr pp)) X (setq f (math-add (math-mul x f) d) X d (math-add (math-mul x d) b) X b (math-add (math-mul x b) (car pp)) X err (math-add (math-abs-approx b) (math-mul abx err)))) X (math-lessp (calcFunc-scf err (- -2 calc-internal-prec)) X (math-abs-approx b))) X (or (not (math-zerop d)) X (not (math-zerop f)) X (progn X (setq x (math-pow (math-neg b) (list 'frac 1 m))) X nil)) X (let* ((g (math-div d b)) X (g2 (math-sqr g)) X (h (math-sub g2 (math-mul 2 (math-div f b)))) X (sq (math-sqrt X (math-mul (1- m) (math-sub (math-mul m h) g2)))) X (gp (math-add g sq)) X (gm (math-sub g sq))) X (if (math-lessp (calcFunc-abssqr gp) (calcFunc-abssqr gm)) X (setq gp gm)) X (setq dx (math-div m gp) X x1 (math-sub x dx)) X (if (and try-newt X (math-lessp (math-abs-approx dx) X (calcFunc-scf (math-abs-approx x) -3))) X (let ((newt (math-poly-newton-root p x1 7))) X (setq tried-newt t X try-newt nil) X (if (math-zerop (cdr newt)) X (setq x (car newt) x1 x) X (if (math-lessp (cdr newt) '(float 1 -6)) X (let ((newt2 (math-poly-newton-root X p (car newt) 20))) X (if (math-zerop (cdr newt2)) X (setq x (car newt2) x1 x) X (setq x (car newt)))))))) X (not (or (eq x x1) X (math-nearly-equal x x1)))) X (let ((cdx (math-abs-approx dx))) X (setq x x1 X tried-newt nil) X (prog1 X (or (<= iters 6) X (math-lessp cdx dxold) X (progn X (if polish X (let ((digs (calcFunc-xpon X (math-div (math-abs-approx x) cdx)))) X (calc-record-why X "*Could not attain full precision") X (if (natnump digs) X (let ((calc-internal-prec (max 3 digs))) X (setq x (math-normalize x)))))) X nil)) X (setq dxold cdx))) X (or polish X (math-lessp (calcFunc-scf (math-abs-approx x) X (- calc-internal-prec)) X dxold)))) X (or (and (math-floatp x) X (math-poly-integer-root x)) X x)) ) X (defun math-solve-above-dummy (x) X (and (not (Math-primp x)) X (if (and (equal (nth 1 x) '(var SOLVEDUM SOLVEDUM)) X (= (length x) 2)) X x X (let ((res nil)) X (while (and (setq x (cdr x)) X (not (setq res (math-solve-above-dummy (car x)))))) X res))) ) X (defun math-solve-find-root-term (x neg) ; sets "t2", "t3" X (if (math-solve-find-root-in-prod x) X (setq t3 neg X t1 x) X (and (memq (car-safe x) '(+ -)) X (or (math-solve-find-root-term (nth 1 x) neg) X (math-solve-find-root-term (nth 2 x) X (if (eq (car x) '-) (not neg) neg))))) ) X (defun math-solve-find-root-in-prod (x) X (and (consp x) X (math-expr-contains x solve-var) X (or (and (eq (car x) 'calcFunc-sqrt) X (setq t2 2)) X (and (eq (car x) '^) X (or (and (memq (math-quarter-integer (nth 2 x)) '(1 2 3)) X (setq t2 2)) X (and (eq (car-safe (nth 2 x)) 'frac) X (eq (nth 2 (nth 2 x)) 3) X (setq t2 3)))) X (and (memq (car x) '(* /)) X (or (and (not (math-expr-contains (nth 1 x) solve-var)) X (math-solve-find-root-in-prod (nth 2 x))) X (and (not (math-expr-contains (nth 2 x) solve-var)) X (math-solve-find-root-in-prod (nth 1 x))))))) ) X X (defun math-solve-system (exprs solve-vars solve-full) X (setq exprs (mapcar 'list (if (Math-vectorp exprs) X (cdr exprs) X (list exprs))) X solve-vars (if (Math-vectorp solve-vars) X (cdr solve-vars) X (list solve-vars))) X (or (let ((math-solve-simplifying nil)) X (math-solve-system-rec exprs solve-vars nil)) X (let ((math-solve-simplifying t)) X (math-solve-system-rec exprs solve-vars nil))) ) X ;;; The following backtracking solver works by choosing a variable ;;; and equation, and trying to solve the equation for the variable. ;;; If it succeeds it calls itself recursively with that variable and ;;; equation removed from their respective lists, and with the solution ;;; added to solns as well as being substituted into all existing ;;; equations. The algorithm terminates when any solution path ;;; manages to remove all the variables from var-list. X ;;; To support calcFunc-roots, entries in eqn-list and solns are ;;; actually lists of equations. X (defun math-solve-system-rec (eqn-list var-list solns) X (if var-list X (let ((v var-list) X (res nil)) X X ;; Try each variable in turn. X (while X (and X v X (let* ((vv (car v)) X (e eqn-list) X (elim (eq (car-safe vv) 'calcFunc-elim))) X (if elim X (setq vv (nth 1 vv))) X X ;; Try each equation in turn. X (while X (and X e X (let ((e2 (car e)) X (eprev nil) X res2) X (setq res nil) X X ;; Try to solve for vv the list of equations e2. X (while (and e2 X (setq res2 (or (and (eq (car e2) eprev) X res2) X (math-solve-for (car e2) 0 vv X solve-full)))) X (setq eprev (car e2) X res (cons (if (eq solve-full 'all) X (cdr res2) X (list res2)) X res) X e2 (cdr e2))) X (if e2 X (setq res nil) X X ;; Found a solution. Now try other variables. X (setq res (nreverse res) X res (math-solve-system-rec X (mapcar X 'math-solve-system-subst X (delq (car e) X (copy-sequence eqn-list))) X (delq (car v) (copy-sequence var-list)) X (let ((math-solve-simplifying nil) X (s (mapcar X (function X (lambda (x) X (cons X (car x) X (math-solve-system-subst X (cdr x))))) X solns))) X (if elim X s X (cons (cons vv (apply 'append res)) X s))))) X (not res)))) X (setq e (cdr e))) X (not res))) X (setq v (cdr v))) X res) X X ;; Eliminated all variables, so now put solution into the proper format. X (setq solns (sort solns X (function X (lambda (x y) X (not (memq (car x) (memq (car y) solve-vars))))))) X (if (eq solve-full 'all) X (math-transpose X (math-normalize X (cons 'vec X (if solns X (mapcar (function (lambda (x) (cons 'vec (cdr x)))) solns) X (mapcar (function (lambda (x) (cons 'vec x))) eqn-list))))) X (math-normalize X (cons 'vec X (if solns SHAR_EOF true || echo 'restore of calc-alg-2.el failed' fi echo 'End of part 5' echo 'File calc-alg-2.el is continued in part 6' echo 6 > _shar_seq_.tmp exit 0 exit 0 # Just in case... -- Kent Landfield INTERNET: kent@sparky.IMD.Sterling.COM Sterling Software, IMD UUCP: uunet!sparky!kent Phone: (402) 291-8300 FAX: (402) 291-4362 Please send comp.sources.misc-related mail to kent@uunet.uu.net.